home *** CD-ROM | disk | FTP | other *** search
/ SGI Freeware 1998 November / Freeware November 1998.img / dist / fw_emacs.idb / usr / freeware / share / emacs / 19.34 / lisp / float.el.z / float.el
Lisp/Scheme  |  1998-10-27  |  15KB  |  459 lines
  1. ;;; float.el --- floating point arithmetic package.
  2.  
  3. ;; Copyright (C) 1986 Free Software Foundation, Inc.
  4.  
  5. ;; Author: Bill Rosenblatt
  6. ;; Maintainer: FSF
  7. ;; Keywords: extensions
  8.  
  9. ;; This file is part of GNU Emacs.
  10.  
  11. ;; GNU Emacs is free software; you can redistribute it and/or modify
  12. ;; it under the terms of the GNU General Public License as published by
  13. ;; the Free Software Foundation; either version 2, or (at your option)
  14. ;; any later version.
  15.  
  16. ;; GNU Emacs is distributed in the hope that it will be useful,
  17. ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  18. ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19. ;; GNU General Public License for more details.
  20.  
  21. ;; You should have received a copy of the GNU General Public License
  22. ;; along with GNU Emacs; see the file COPYING.  If not, write to the
  23. ;; Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  24. ;; Boston, MA 02111-1307, USA.
  25.  
  26. ;;; Commentary:
  27.  
  28. ;; Floating point numbers are represented by dot-pairs (mant . exp)
  29. ;; where mant is the 24-bit signed integral mantissa and exp is the
  30. ;; base 2 exponent.
  31. ;;
  32. ;; Emacs LISP supports a 24-bit signed integer data type, which has a
  33. ;; range of -(2**23) to +(2**23)-1, or -8388608 to 8388607 decimal.
  34. ;; This gives six significant decimal digit accuracy.  Exponents can
  35. ;; be anything in the range -(2**23) to +(2**23)-1.
  36. ;;
  37. ;; User interface:
  38. ;; function f converts from integer to floating point
  39. ;; function string-to-float converts from string to floating point
  40. ;; function fint converts a floating point to integer (with truncation)
  41. ;; function float-to-string converts from floating point to string
  42. ;;                   
  43. ;; Caveats:
  44. ;; -  Exponents outside of the range of +/-100 or so will cause certain 
  45. ;;    functions (especially conversion routines) to take forever.
  46. ;; -  Very little checking is done for fixed point overflow/underflow.
  47. ;; -  No checking is done for over/underflow of the exponent
  48. ;;    (hardly necessary when exponent can be 2**23).
  49. ;; 
  50. ;;
  51. ;; Bill Rosenblatt
  52. ;; June 20, 1986
  53. ;;
  54.  
  55. ;;; Code:
  56.  
  57. ;; fundamental implementation constants
  58. (defconst exp-base 2
  59.   "Base of exponent in this floating point representation.")
  60.  
  61. (defconst mantissa-bits 24
  62.   "Number of significant bits in this floating point representation.")
  63.  
  64. (defconst decimal-digits 6
  65.   "Number of decimal digits expected to be accurate.")
  66.  
  67. (defconst expt-digits 2
  68.   "Maximum permitted digits in a scientific notation exponent.")
  69.  
  70. ;; other constants
  71. (defconst maxbit (1- mantissa-bits)
  72.   "Number of highest bit")
  73.  
  74. (defconst mantissa-maxval (1- (ash 1 maxbit))
  75.   "Maximum permissible value of mantissa")
  76.  
  77. (defconst mantissa-minval (ash 1 maxbit)
  78.   "Minimum permissible value of mantissa")
  79.  
  80. (defconst floating-point-regexp
  81.   "^[ \t]*\\(-?\\)\\([0-9]*\\)\
  82. \\(\\.\\([0-9]*\\)\\|\\)\
  83. \\(\\(\\([Ee]\\)\\(-?\\)\\([0-9][0-9]*\\)\\)\\|\\)[ \t]*$"
  84.   "Regular expression to match floating point numbers.  Extract matches:
  85. 1 - minus sign
  86. 2 - integer part
  87. 4 - fractional part
  88. 8 - minus sign for power of ten
  89. 9 - power of ten
  90. ")
  91.  
  92. (defconst high-bit-mask (ash 1 maxbit)
  93.   "Masks all bits except the high-order (sign) bit.")
  94.  
  95. (defconst second-bit-mask (ash 1 (1- maxbit))
  96.   "Masks all bits except the highest-order magnitude bit")
  97.  
  98. ;; various useful floating point constants
  99. (setq _f0 '(0 . 1))
  100.  
  101. (setq _f1/2 '(4194304 . -23))
  102.  
  103. (setq _f1 '(4194304 . -22))
  104.  
  105. (setq _f10 '(5242880 . -19))
  106.  
  107. ;; support for decimal conversion routines
  108. (setq powers-of-10 (make-vector (1+ decimal-digits) _f1))
  109. (aset powers-of-10 1 _f10)
  110. (aset powers-of-10 2 '(6553600 . -16))
  111. (aset powers-of-10 3 '(8192000 . -13))
  112. (aset powers-of-10 4 '(5120000 . -9))
  113. (aset powers-of-10 5 '(6400000 . -6))
  114. (aset powers-of-10 6 '(8000000 . -3))
  115.  
  116. (setq all-decimal-digs-minval (aref powers-of-10 (1- decimal-digits))
  117.       highest-power-of-10 (aref powers-of-10 decimal-digits))
  118.  
  119. (defun fashl (fnum)            ; floating-point arithmetic shift left
  120.   (cons (ash (car fnum) 1) (1- (cdr fnum))))
  121.  
  122. (defun fashr (fnum)            ; floating point arithmetic shift right
  123.   (cons (ash (car fnum) -1) (1+ (cdr fnum))))
  124.  
  125. (defun normalize (fnum)
  126.   (if (> (car fnum) 0)            ; make sure next-to-highest bit is set
  127.       (while (zerop (logand (car fnum) second-bit-mask))
  128.     (setq fnum (fashl fnum)))
  129.     (if (< (car fnum) 0)        ; make sure highest bit is set
  130.     (while (zerop (logand (car fnum) high-bit-mask))
  131.       (setq fnum (fashl fnum)))
  132.       (setq fnum _f0)))            ; "standard 0"
  133.   fnum)
  134.       
  135. (defun abs (n)                ; integer absolute value
  136.   (if (>= n 0) n (- n)))
  137.  
  138. (defun fabs (fnum)            ; re-normalize after taking abs value
  139.   (normalize (cons (abs (car fnum)) (cdr fnum))))
  140.  
  141. (defun xor (a b)            ; logical exclusive or
  142.   (and (or a b) (not (and a b))))
  143.  
  144. (defun same-sign (a b)            ; two f-p numbers have same sign?
  145.   (not (xor (natnump (car a)) (natnump (car b)))))
  146.  
  147. (defun extract-match (str i)        ; used after string-match
  148.   (condition-case ()
  149.       (substring str (match-beginning i) (match-end i))
  150.     (error "")))
  151.  
  152. ;; support for the multiplication function
  153. (setq halfword-bits (/ mantissa-bits 2)    ; bits in a halfword
  154.       masklo (1- (ash 1 halfword-bits)) ; isolate the lower halfword
  155.       maskhi (lognot masklo)        ; isolate the upper halfword
  156.       round-limit (ash 1 (/ halfword-bits 2)))
  157.  
  158. (defun hihalf (n)            ; return high halfword, shifted down
  159.   (ash (logand n maskhi) (- halfword-bits)))
  160.  
  161. (defun lohalf (n)            ; return low halfword
  162.   (logand n masklo))
  163.  
  164. ;; Visible functions
  165.  
  166. ;; Arithmetic functions
  167. (defun f+ (a1 a2)
  168.   "Returns the sum of two floating point numbers."
  169.   (let ((f1 (fmax a1 a2))
  170.     (f2 (fmin a1 a2)))
  171.     (if (same-sign a1 a2)
  172.     (setq f1 (fashr f1)        ; shift right to avoid overflow
  173.           f2 (fashr f2)))
  174.     (normalize
  175.      (cons (+ (car f1) (ash (car f2) (- (cdr f2) (cdr f1))))
  176.        (cdr f1)))))
  177.  
  178. (defun f- (a1 &optional a2)        ; unary or binary minus
  179.   "Returns the difference of two floating point numbers."
  180.   (if a2
  181.       (f+ a1 (f- a2))
  182.     (normalize (cons (- (car a1)) (cdr a1)))))
  183.  
  184. (defun f* (a1 a2)            ; multiply in halfword chunks
  185.   "Returns the product of two floating point numbers."
  186.   (let* ((i1 (car (fabs a1)))
  187.      (i2 (car (fabs a2)))
  188.      (sign (not (same-sign a1 a2)))
  189.      (prodlo (+ (hihalf (* (lohalf i1) (lohalf i2)))
  190.             (lohalf (* (hihalf i1) (lohalf i2)))
  191.             (lohalf (* (lohalf i1) (hihalf i2)))))
  192.      (prodhi (+ (* (hihalf i1) (hihalf i2))
  193.             (hihalf (* (hihalf i1) (lohalf i2)))
  194.             (hihalf (* (lohalf i1) (hihalf i2)))
  195.             (hihalf prodlo))))
  196.     (if (> (lohalf prodlo) round-limit)
  197.     (setq prodhi (1+ prodhi)))    ; round off truncated bits
  198.     (normalize
  199.      (cons (if sign (- prodhi) prodhi)
  200.        (+ (cdr (fabs a1)) (cdr (fabs a2)) mantissa-bits)))))
  201.  
  202. (defun f/ (a1 a2)            ; SLOW subtract-and-shift algorithm
  203.   "Returns the quotient of two floating point numbers."
  204.   (if (zerop (car a2))            ; if divide by 0
  205.       (signal 'arith-error (list "attempt to divide by zero" a1 a2))
  206.     (let ((bits (1- maxbit))
  207.       (quotient 0) 
  208.       (dividend (car (fabs a1)))
  209.       (divisor (car (fabs a2)))
  210.       (sign (not (same-sign a1 a2))))
  211.       (while (natnump bits)
  212.     (if (< (- dividend divisor) 0)
  213.         (setq quotient (ash quotient 1))
  214.       (setq quotient (1+ (ash quotient 1))
  215.         dividend (- dividend divisor)))
  216.     (setq dividend (ash dividend 1)
  217.           bits (1- bits)))
  218.       (normalize
  219.        (cons (if sign (- quotient) quotient)
  220.          (- (cdr (fabs a1)) (cdr (fabs a2)) (1- maxbit)))))))
  221.   
  222. (defun f% (a1 a2)
  223.   "Returns the remainder of first floating point number divided by second."
  224.   (f- a1 (f* (ftrunc (f/ a1 a2)) a2)))
  225.       
  226.  
  227. ;; Comparison functions
  228. (defun f= (a1 a2)
  229.   "Returns t if two floating point numbers are equal, nil otherwise."
  230.   (equal a1 a2))
  231.  
  232. (defun f> (a1 a2)
  233.   "Returns t if first floating point number is greater than second,
  234. nil otherwise."
  235.   (cond ((and (natnump (car a1)) (< (car a2) 0)) 
  236.      t)                ; a1 nonnegative, a2 negative
  237.     ((and (> (car a1) 0) (<= (car a2) 0))
  238.      t)                ; a1 positive, a2 nonpositive
  239.     ((and (<= (car a1) 0) (natnump (car a2)))
  240.      nil)                ; a1 nonpos, a2 nonneg
  241.     ((/= (cdr a1) (cdr a2))        ; same signs.  exponents differ
  242.      (> (cdr a1) (cdr a2)))        ; compare the mantissas.
  243.     (t
  244.      (> (car a1) (car a2)))))    ; same exponents.
  245.  
  246. (defun f>= (a1 a2)
  247.   "Returns t if first floating point number is greater than or equal to 
  248. second, nil otherwise."
  249.   (or (f> a1 a2) (f= a1 a2)))
  250.  
  251. (defun f< (a1 a2)
  252.   "Returns t if first floating point number is less than second,
  253. nil otherwise."
  254.   (not (f>= a1 a2)))
  255.  
  256. (defun f<= (a1 a2)
  257.   "Returns t if first floating point number is less than or equal to
  258. second, nil otherwise."
  259.   (not (f> a1 a2)))
  260.  
  261. (defun f/= (a1 a2)
  262.   "Returns t if first floating point number is not equal to second,
  263. nil otherwise."
  264.   (not (f= a1 a2)))
  265.  
  266. (defun fmin (a1 a2)
  267.   "Returns the minimum of two floating point numbers."
  268.   (if (f< a1 a2) a1 a2))
  269.  
  270. (defun fmax (a1 a2)
  271.   "Returns the maximum of two floating point numbers."
  272.   (if (f> a1 a2) a1 a2))
  273.       
  274. (defun fzerop (fnum)
  275.   "Returns t if the floating point number is zero, nil otherwise."
  276.   (= (car fnum) 0))
  277.  
  278. (defun floatp (fnum)
  279.   "Returns t if the arg is a floating point number, nil otherwise."
  280.   (and (consp fnum) (integerp (car fnum)) (integerp (cdr fnum))))
  281.  
  282. ;; Conversion routines
  283. (defun f (int)
  284.   "Convert the integer argument to floating point, like a C cast operator."
  285.   (normalize (cons int '0)))
  286.  
  287. (defun int-to-hex-string (int)
  288.   "Convert the integer argument to a C-style hexadecimal string."
  289.   (let ((shiftval -20)
  290.     (str "0x")
  291.     (hex-chars "0123456789ABCDEF"))
  292.     (while (<= shiftval 0)
  293.       (setq str (concat str (char-to-string 
  294.             (aref hex-chars
  295.                   (logand (lsh int shiftval) 15))))
  296.         shiftval (+ shiftval 4)))
  297.     str))
  298.  
  299. (defun ftrunc (fnum)            ; truncate fractional part
  300.   "Truncate the fractional part of a floating point number."
  301.   (cond ((natnump (cdr fnum))        ; it's all integer, return number as is
  302.      fnum)
  303.     ((<= (cdr fnum) (- maxbit))    ; it's all fractional, return 0
  304.      '(0 . 1))
  305.     (t                ; otherwise mask out fractional bits
  306.      (let ((mant (car fnum)) (exp (cdr fnum)))
  307.        (normalize 
  308.         (cons (if (natnump mant)    ; if negative, use absolute value
  309.               (ash (ash mant exp) (- exp))
  310.             (- (ash (ash (- mant) exp) (- exp))))
  311.           exp))))))
  312.  
  313. (defun fint (fnum)            ; truncate and convert to integer
  314.   "Convert the floating point number to integer, with truncation, 
  315. like a C cast operator."
  316.   (let* ((tf (ftrunc fnum)) (tint (car tf)) (texp (cdr tf)))
  317.     (cond ((>= texp mantissa-bits)    ; too high, return "maxint"
  318.        mantissa-maxval)
  319.       ((<= texp (- mantissa-bits))    ; too low, return "minint"
  320.        mantissa-minval)
  321.       (t                ; in range
  322.        (ash tint texp)))))        ; shift so that exponent is 0
  323.  
  324. (defun float-to-string (fnum &optional sci)
  325.   "Convert the floating point number to a decimal string.
  326. Optional second argument non-nil means use scientific notation."
  327.   (let* ((value (fabs fnum)) (sign (< (car fnum) 0))
  328.      (power 0) (result 0) (str "") 
  329.      (temp 0) (pow10 _f1))
  330.  
  331.     (if (f= fnum _f0)
  332.     "0"
  333.       (if (f>= value _f1)            ; find largest power of 10 <= value
  334.       (progn                ; value >= 1, power is positive
  335.         (while (f<= (setq temp (f* pow10 highest-power-of-10)) value)
  336.           (setq pow10 temp
  337.             power (+ power decimal-digits)))
  338.         (while (f<= (setq temp (f* pow10 _f10)) value)
  339.           (setq pow10 temp
  340.             power (1+ power))))
  341.     (progn                ; value < 1, power is negative
  342.       (while (f> (setq temp (f/ pow10 highest-power-of-10)) value)
  343.         (setq pow10 temp
  344.           power (- power decimal-digits)))
  345.       (while (f> pow10 value)
  346.         (setq pow10 (f/ pow10 _f10)
  347.           power (1- power)))))
  348.                       ; get value in range 100000 to 999999
  349.       (setq value (f* (f/ value pow10) all-decimal-digs-minval)
  350.         result (ftrunc value))
  351.       (let (int)
  352.     (if (f> (f- value result) _f1/2)    ; round up if remainder > 0.5
  353.         (setq int (1+ (fint result)))
  354.       (setq int (fint result)))
  355.     (setq str (int-to-string int))
  356.     (if (>= int 1000000)
  357.         (setq power (1+ power))))
  358.  
  359.       (if sci                ; scientific notation
  360.       (setq str (concat (substring str 0 1) "." (substring str 1)
  361.                 "E" (int-to-string power)))
  362.  
  363.                       ; regular decimal string
  364.     (cond ((>= power (1- decimal-digits))
  365.                       ; large power, append zeroes
  366.            (let ((zeroes (- power decimal-digits)))
  367.          (while (natnump zeroes)
  368.            (setq str (concat str "0")
  369.              zeroes (1- zeroes)))))
  370.  
  371.                       ; negative power, prepend decimal
  372.           ((< power 0)        ; point and zeroes
  373.            (let ((zeroes (- (- power) 2)))
  374.          (while (natnump zeroes)
  375.            (setq str (concat "0" str)
  376.              zeroes (1- zeroes)))
  377.          (setq str (concat "0." str))))
  378.  
  379.           (t                ; in range, insert decimal point
  380.            (setq str (concat
  381.               (substring str 0 (1+ power))
  382.               "."
  383.               (substring str (1+ power)))))))
  384.  
  385.       (if sign                ; if negative, prepend minus sign
  386.       (concat "-" str)
  387.     str))))
  388.  
  389.     
  390. ;; string to float conversion.
  391. ;; accepts scientific notation, but ignores anything after the first two
  392. ;; digits of the exponent.
  393. (defun string-to-float (str)
  394.   "Convert the string to a floating point number.
  395. Accepts a decimal string in scientific notation,  with exponent preceded
  396. by either E or e.  Only the six most significant digits of the integer
  397. and fractional parts are used; only the first two digits of the exponent
  398. are used.  Negative signs preceding both the decimal number and the exponent
  399. are recognized."
  400.  
  401.   (if (string-match floating-point-regexp str 0)
  402.       (let (power)
  403.     (f*
  404.      ; calculate the mantissa
  405.      (let* ((int-subst (extract-match str 2))
  406.         (fract-subst (extract-match str 4))
  407.         (digit-string (concat int-subst fract-subst))
  408.         (mant-sign (equal (extract-match str 1) "-"))
  409.         (leading-0s 0) (round-up nil))
  410.  
  411.        ; get rid of leading 0's
  412.        (setq power (- (length int-subst) decimal-digits))
  413.        (while (and (< leading-0s (length digit-string))
  414.                (= (aref digit-string leading-0s) ?0))
  415.          (setq leading-0s (1+ leading-0s)))
  416.        (setq power (- power leading-0s)
  417.          digit-string (substring digit-string leading-0s))
  418.        
  419.        ; if more than 6 digits, round off
  420.        (if (> (length digit-string) decimal-digits)
  421.            (setq round-up (>= (aref digit-string decimal-digits) ?5)
  422.              digit-string (substring digit-string 0 decimal-digits))
  423.          (setq power (+ power (- decimal-digits (length digit-string)))))
  424.  
  425.        ; round up and add minus sign, if necessary
  426.        (f (* (+ (string-to-int digit-string)
  427.             (if round-up 1 0))
  428.          (if mant-sign -1 1))))
  429.        
  430.      ; calculate the exponent (power of ten)
  431.      (let* ((expt-subst (extract-match str 9))
  432.         (expt-sign (equal (extract-match str 8) "-"))
  433.         (expt 0) (chunks 0) (tens 0) (exponent _f1)
  434.         (func 'f*))
  435.  
  436.        (setq expt (+ (* (string-to-int
  437.                  (substring expt-subst 0
  438.                     (min expt-digits (length expt-subst))))
  439.                 (if expt-sign -1 1))
  440.              power))
  441.        (if (< expt 0)        ; if power of 10 negative
  442.            (setq expt (- expt)    ; take abs val of exponent
  443.              func 'f/))        ; and set up to divide, not multiply
  444.  
  445.        (setq chunks (/ expt decimal-digits)
  446.          tens (% expt decimal-digits))
  447.        ; divide or multiply by "chunks" of 10**6
  448.        (while (> chunks 0)    
  449.          (setq exponent (funcall func exponent highest-power-of-10)
  450.            chunks (1- chunks)))
  451.        ; divide or multiply by remaining power of ten
  452.        (funcall func exponent (aref powers-of-10 tens)))))
  453.           
  454.     _f0))                ; if invalid, return 0
  455.  
  456. (provide 'float)
  457.  
  458. ;;; float.el ends here
  459.